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1. Introduction 

The solar magnetic field is an important quantity which couples the so- 
lar interior with the photosphere and atmosphere. Knowledge regarding 
the coronal magnetic field plays a key role for eruptive phenomena, e.g. 
coronal mass ejection, flares and eruptive prominences. Unfortunately a 
direct measurement of the coronal magnetic field is extremely difficult. 
In principle one can use the polarization of emissions from magnetic 
sensitive coronal line transitions to draw conclusions about the coronal 
magnetic field. These lines, however, are very faint so that in the past 
they have only occasionally been observed (e.g., (House, 1977; Arnaud 
and Newkirk, 1987; Judge, 1998)). In a recent study (Judge et al., 
2001) conclude that several forbidden lines (e.g. in Fe XIII, He I, Mg 
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VIII and Si IX) may be used to determine the coronal magnetic field. 
They further concluded that space born missions are not needed for 
such kind of coronal magnetometers but a high, dry mountain site. In 
their study they propose a focal plane instrument devoted to the l/im 
region. These authors also point out that besides the observational 
part, a further major problem is the interpretation of the data. The 
line-of-sight integration inherent in these observations make the data 
analysis a badly posed inversion problem. Presently, algorithms based 
on vector tomography are studied to find out to which extend the 
coronal magnetic field can be reconstructed from these observations 
(Maxim Kramar and Bernd Inhester, private communication). 

Despite these promising new developments regarding the principle 
possibility of coronal B-field measurements we have to face the fact that 
currently and probably also in near future high quality direct measure- 
ments of the coronal magnetic field are not available. As an alternative 
to the above measurements, a number of authors have modelled the 
coronal magnetic field by extrapolation from more sound photospheric 
magnetic field observations. 

It is generally assumed that the magnetic pressure in the corona 
is much higher than the plasma pressure (small plasma 0) and that 
therefore the magnetic field is nearly force- free (for a critical view of this 
assumption see (Gary, 2001)). The extrapolation methods based on this 
assumption include potential field extrapolation (Schmidt, 1964; Semel, 
1967), linear force- free field extrapolation (Chiu and Hilton, 1977; See- 
hafer, 1978; Seehafer, 1982; Semel, 1988) and nonlinear force-free field 
extrapolation (Amari et al., 1997). Methods for the extrapolation non 
force-free fields and have been developed by (Petrie and Neukirch, 
2000; Wiegelmann and Inhester, 2003). Potential fields can be de- 
termined directly from line-of-sight magnetogram data (e.g. MDI on 
SOHO). Linear force-free fields can as well be calculated from line-of- 
sight magnetograms, but contain a free parameter a which has to be 
computed from additional data, e.g. with fitting procedures which try 
to match the field model with observed coronal plasma loops using data 
from e.g. EIT (Wiegelmann and Neukirch, 2002) or Yohkoh (Carcedo 
et al., 2003). Unfortunately potential fields and linear force- free fields 
do not contain free energy and are very probably a poor approximation 
for an active region prior to an eruption. By free energy we understand 
energy which can be released during an eruption. A linear force free 
field has more energy than a potential field. This energy can, however, 
not be released during an eruption related to ideal or resistive MHD- 
instabilities because a linear force free field cannot rapidly relax to 
a potential field. The reason is that the magnetic helicity is strictly 
conserved for ideal MHD and approximately conserved for resistive 
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processes. (The magnetic helicity is dissipated slower than the magnetic 
energy, see (Berger, 1984)). A nonlinear force free field can, however, 
relax to a linear force free field with the same magnetic helicity. In 
this sense a nonlinear force free field has free energy available for an 
eruption. Consequently investigations regarding non-linear force-free 
fields are essential to understand eruptive phenomena. 

The calculation of non linear force-free fields is complicated by the 
intrinsic nonlinearity of the underlying mathematical problem. From 
the observational point of view the non linear reconstruction is also 
more challenging because photospheric vector magnetograph data are 
required. Unfortunately the transversal component of the photospheric 
B-field is measured with significant lower accuracy as the line of sight 
component. An additional problem is that the transversal magnetic 
field is only known with respect to an 180° ambiguity and a prepro- 
cessing of the raw data is necessary to resolve this ambiguity (One 
possibility is the minimum energy method by (Metcalf, 1994) used for 
vector magnetogram data from IVM in Hawaii.). 

Several methods have been proposed to compute nonlinear force- free 
fields: 

— A conceptionally simple method is to reformulate the force-free 
equations (2-3) in such way that they can be used for an upward 
integration of the vector magnetogram into the corona (Wu et 
al., 1990; Amari et al., 1997). This direct extrapolation is an ill 
posed problem for the elliptic equations (2-3) and consequently the 
method is limited to low heights. In particular one finds that an 
erroneous exponential growth of the magnetic field with increasing 
height is a typical behaviour. 

— An alternative approach is to use a Grad- Rubin method ((Sakurai, 
1981; Amari et al.1999)). This method uses a potential field as 
initial equilibrium and then progressively currents are introduced 
into the system and the fields are relaxed towards a force-free state. 
The method is especial useful for small deviations from a poten- 
tial field with small values of a and modest non linearities. The 
method requires an explicit calculation of the a distribution on the 
photosphere. In principle the computation of a is straight forward 
a(x,y) = dBy ^ d ^\xy) ^ dV this is inaccurate for observational 
data for the following reasons. First one needs the transversal 
component of the photospheric magnetic field which is measured 
with lower accuracy than the line of sight magnetic field. In a 
second step one has to take the horizontal derivatives (x,y)of these 
inaccurate values and finally one has to divide through the normal 
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magnetic field B z which causes additional problems where \B Z \ is 
small. These errors cumulate in the photospheric a distribution. 



— A third possibility is to use the method of MHD relaxation (Chodura 
and Schliiter, 1981; Roumeliotis, 1996). The idea is to start with a 
suitable magnetic field which is not in equilibrium and to relax it 
into a force-free state. For test configurations (Low and Lou, 1989) 
the MHD relaxation method converges to the exact solution, but 
with less accuracy than the optimization approach discussed below 
(Wiegelmann and Neukirch, 2003). 

— In the optimization approach (Wheatland et al., 2000) a functional 
containing the force-free equations is minimized. The method di- 
rectly uses the measured vector magnetograph data and an ex- 
plicit computation of a is not necessary. Another advantage of 
the method is that the quality of the reconstructed magnetic field 
(force-free and solenoidal condition) is controlled automatically 
within the iteration procedure. A difficulty of the method is that it 
requires boundary conditions on all boundaries of a computational 
box while for observational data only the bottom boundary data 
are known. Within this paper we are dealing with this problem 
and extend the optimization method accordingly. 

While the low (3 plasma in the lower corona can be described with 
the non linear force-free approach, it is necessary to include plasma 
pressure and solar gravitation to describe regions with a finite plasma (3, 
e.g. helmet streamers. (Wiegelmann and Inhester, 2003) extended the 
force-free optimization with the aim to include these forces and showed 
that the method converges for test configurations. The difficulties re- 
garding the lateral and top boundaries are analog to the non-linear 
force-free case. The non force-free reconstruction requires additional 
data as input, e.g. the coronal plasma density distribution from an 
assumed model or computed with help of tomographic methods. 

We outline the paper as follows. In section 2 we provide the basic 
equations of the modified optimization method and derive the iteration 
equations. In section 3 we specify useful forms of the weighting function 
in the boundary regions. Section 4 contains test runs regarding the non 
linear force-free case and section 5 consistency checks for non force-free 
configurations. We draw conclusions in section 6 and give an outlook 
for further research. 
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2. Basic equations 

Force-free coronal magnetic fields have to obey the equations 

jxB = 0, (1) 
V x B = fi j, (2) 
V-B = 0. (3) 

The force-free approach is valid in the low corona where the plasma (3 
is small. For extended structures, e.g. helmet streamers the plasma (3 
increases and the force-free assumption is not valid anymore. Therefore 
it is necessary to consider the effect of plasma pressure and gravity here 
and solve the magneto hydro static equations (MHS). 

j x B - VP - pW = 0, (4) 
V x B = p j, (5) 
V-B = 0, (6) 

where B is the magnetic field, j the electric current density, P the 
plasma pressure, p the plasma density, po the vacuum permeability 
and ^ the solar gravity potential. We define the functional 

L = [ w(x,y,z) B 2 (n 2 a + n 2 b )d 3 x, (7) 
Jv 

with 



B- 2 [(V x B)] 
B- 2 [(VxB)xB- /i (VP + pV*)] (MHS) 



(force-free fields) 



Q h = B- 2 [(V-B)B]. 



(8) 
(9) 



w(x, y, z) is a weighting function. Useful forms of the weighting function 
will be discussed below. 

For the force-free case the functional is given explicitly as 



v 



w(x,y,z) B~ 2 |(V x B) x B| 2 + |V • B| 2 d 3 x, 



(10) 



and it is obvious that (for w > 0) the force-free equations (1-3) are 
fulfilled when L is equal zero. For the non force- free case L is given 
explicitly as 



L 



= f w( 
Jv 



x,y,z) 



B _2 |(VxB) xB- M (VP + pV^)| 2 + |V-B| 



cfx, 



(11) 
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and when the functional reaches (for w > 0) its minimum at L = 
then the MHS equations (4-6) are fulfilled. 

The following discussion is equivalent for the force-free and non 
force-free case. Without weighting function (w = 1) the method has 
been developed by (Wheatland et al., 2000) for the force- free case and 
by (Wiegelmann and Inhester, 2003) for the non force-free case. For 
w(x,y,z) = 1 the optimization method requires that the magnetic field 
is given on all (6 for a rectangular computational box) boundaries. This 
causes a serious limitation of the method because such data are only 
available for model configurations. For the reconstruction of the coronal 
magnetic field it is necessary to develop a method which reconstructs 
the magnetic field only from photospheric vector magnetograms. Vec- 
tor magnetograms provide boundary conditions only for the bottom 
boundary of a computational box while the other five boundaries re- 
main unknown. Without a weighting function all six boundaries of 
the computational box have equal rights and influence the solution 
in the box. It is therefore important to diminish the effect of the top 
and lateral boundaries on the magnetic field inside the computational 
box. This can be done either by including a variation of B not only 
in the interior but also on those boundaries where B is unknown. This 
approach, however, is numerically difficult because it involves two types 
of variations. We show that it is essentially equivalent to introducing 
finite size boundary regions on those boundaries where B is unknown 
with the weighting function w(x,y,z) different from unity. 

The idea is do define an interior physical region where we want to 
calculate the magnetic field so that it fulfills the force-free or MHS equa- 
tions. This region is in the center of the box (including the photosphere) 
with w = 1. The computational box additionally includes boundary 
layers towards the lateral and top boundary where w decreases to 
at the computational boundary. Consequently the method weights 
deviations from the force-free state (or MHS-state) less severely close 
to the boundary. The use of a weighting function has been proposed for 
the force- free case by (Wheatland et al., 2000) in the conclusions but 
no iteration equations or test simulations have been presented. Here 
we provide these iteration equations for the more generalized case. 
We carry out several tests to investigate the optimum shape of the 
weighting function and how the size of the boundary layer influences 
the quality of the reconstruction. 

We minimize equation (7) with respect to an iteration parameter t 
(see Appendix A for details) and obtain an iteration equation for the 
magnetic field 

/»Ffc-/».GA (12) 
2 dt Jv dt Js dt y ' 
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F = w F + (fi a x B) x Vw + (fi b • B) Vw 
G = wG 



(13) 
(14) 



F = Vx(fi a xB)-n a x(VxB) 

+V(f2 b • B) - ft b (V • B) + {ill + fi b) B 



(15) 



G = n x (O a x B) - n(fi b • B) 



(16) 



and h is the inward unit vector on the surface S. The surface integral 
in (12) vanishes if the magnetic field is described on the boundaries 
of a computational box. Inside the computational box we iterate the 
magnetic field with 



which insures that L is monotonically decreasing. 
2.1. Algorithm 

We compute the 3D-coronal magnetic field in a numerical box using 
the following steps. 

— As a start configuration we use the measured normal component 
B z of the magnetic field to calculate a potential magnetic field 
in the whole box with help of a Fourier representation (Seehafer, 



— For non force-free (finite (3) configurations the plasma density 
distribution is described in the box. This step is unnecessary for 
force-free (/3 <C 1) fields. 

— We use vector magnetograph data to describe the bottom bound- 
ary (photosphere) of the computational box. On the lateral and 
top boundaries the field is chosen from the potential field above. 

— We iterate for the magnetic field inside the computational box 
with (17) using a Landweber-iteration (see e.g. (Louis, 1989)). The 
continuous form of (17) guaranties a monotonically decreasing L. 
This is as well ensured in the discretized form if the iteration step 
dt is sufficiently small. The code checks if L(t + dt) < L{t) after 
each time step. If the condition is not fulfilled, the iteration step 
is repeated with dt reduced by a factor of 2. After each successful 
iteration step we increase dt slowly by a factor of 1.01 to allow 
the time step to become as large as possible with respect to the 




(17) 



1978). 
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stability condition. The iteration stops if dt falls below a limiting 
value, e.g. 1/100 of the initial iteration step 1 in the current version 
of the code. 

Let us remark that the main numerics of the optimization code is 
similar for the method with and without a weighting function. The 
main problem for the optimization method without weighting function 
is that it requires the vector magnetic field on all six boundaries of 
a computational box. As only the bottom boundary is measured one 
has to make assumptions regarding the lateral and top boundary, e.g. 
assume a potential field. In general this leads to inconsistent boundary 
conditions (See (Aly, 1989) regarding the compatibility of photospheric 
vector magnetograph data.) and consequently a bad quality of the 
reconstructed magnetic field. With help of the weighting function the 
five inconsistent boundaries are replaced by boundary layers and con- 
sequently we get more flexible boundaries around the physical domain 
which will adjust automatically during the iteration. The idea of in- 
troducing a boundary layer with w < 1 is to reduce the dependence 
of the solution in the interior of the box from the unknown boundary 
conditions. Since we have no measurements on these boundaries any 
choice of the boundary conditions is a mere guess. The aim is only to 
allow the solution in the interior to evolve more independently from the 
boundary conditions chosen. So the advantage of the boundary layer is 
a higher degree of independence of the solution in the interior from the 
chosen boundary. The price we have to pay is a higher computing time, 
as the magnetic field has to be iterated within the whole computational 
box which includes the physical domain as well as the boundary layers. 



3. Special forms of the weighting function w(x,y,z). 

We want to use the weighting function to deal with the unknown top 
and lateral boundaries. We define an inner physical domain V{ with 
w = 1 and boundary layers Vb where w decreases monotonically from 1 
to through the outer numerical boundary layer with the thickness d. 
Consequently w becomes one-dimensional in each boundary layer (e.g. 
w = w(z) at the top boundary layer) and we get Vt« = n ^. The 
surface integrals vanish on all boundaries because w = on the top 
and lateral boundaries and ^ = on the bottom boundary where the 
magnetic field is measured with vectormagnetographs. Consequently 

1 We find that the time step dt keeps on decreasing recurrently when the solution 
has converged. We never found a further improvement of L after dt has once fallen 
below 1/100 of the initial iteration step. 
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(12) reduces to 

1 dL r dB _ - r dB 



1 dL f dB „ /■ aB ~ o 
- — = - —— • F d 3 x — / 18 

2 9t Jv b dt v ; 

1 dL r dB , 3 /■ 9B 3 
F d 6 x- w — ■ F d d x 



2 dt 9t Jy 6 at 

[(n a x B) x Vw + (« b • B) Vw] d 3 x (19) 



dB 

v b dt 

1 dL f dB „ ,3 / 5B 

2 di 



= - / — • F d 3 x - / «; — • F d 3 x 

JVi dt Jv b dt 

- / ^ ^ • x (Sl a x B) - n (fi b • B)] d 3 x. (20) 
Jv b dn ot 

It is interesting to investigate the limit of an infinitesimally thin bound- 
ary layer d — > in (20). The thinner the boundary layer becomes 
the steeper is ^ and for an infinitesimally thin boundary layer the 
gradient becomes infinity. The boundary layer is constructed in such 
way that independent from the sheet thickness d we have Jq ^ dn = 1 
which remains true also for d —> 0. In the limit of d —> the term 
h x (fi a x B) — h (fib • B) remains constant through the sheet and 
the integration regarding dn can be carried out explicitly. Consequently 
only a surface integral remains as the last term in equation (20). The 
second integral in (20) vanishes for d — ► and we get 

This exactly coincides with the non weighted case. Consequently equa- 
tion (20) is a generalization of the usual optimization equation (21). 



4. Tests for non linear force-free configurations. 

To test our code we use a semi-analytic model active region developed 
by (Low and Lou, 1989). We use the Low and Lou solution with I = 0.5 
and 4> = 1.4 as test. The normal photospheric magnetic field is normal- 
ized to a maximum of 800 Gauss. Figure 2 shows the normal magnetic 
field B z on the photosphere for this configuration. The framed region 
contains the physical domain. We investigate two cases, in Force-Free I 
the physical domain is approximately flux-balanced and in Force-Free 
II not. We are interested in reconstructing an inner region (physical 
domain) of 40 x 40 x 20 points and diagnose Lj and |J x B| aver- 
aged over the physical domain here. We also diagnose L in the whole 
computational box. 
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Figure 1. Different profiles for the weighting function w. The solid lines shows a 
linear profile, the dotted line a quadratic, the dashed line a cos and the dashdotted 
line a tanh profile. All profiles are equal one at the physical boundary, decrease 
monotonically within the boundary layer and reach zero at the boundary of the 
computational box. 




Figure 2. The pictures show artificial magnetograms extracted from the Low and 
Lou solution with I — 0.5 and <f> = 1.4. The framed region corresponds to to the 
physical domain with a resolution of 40 x 40 pixels. The computational box includes 
boundary layers of nd = 10 points towards each boundary. The colour-coding shows 
the normal magnetic field strength on the photosphere. In the left hand picture 
(Force-Free I) the active region is centered and in the right hand picture displaced 
(Force- Free II). The latter simulates data, where a significant amount of flux is not 
balanced. 
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Figure 3. Top row: The left hand side shows some field lines for the original Low and 
Lou configuration (I — 0.5 and (f> = 1.4). The colour-coding shows the normal mag- 
netic field on the photosphere. We choose identical inner footpoints positions for all 
panels. The right hand side contains a corresponding potential field reconstruction. 
The yellow dotted field lines correspond to the original Low and Lou solution in all 
pictures. Center row: The left hand side shows a reconstruction without weighting 
function and the right hand side with a linear weighting function and a 5 pixel 
boundary layer. Bottom row: Both pictures show a reconstruction with cos-profile. 
In the left hand picture a boundary layer of 10 points was used and in the right 
hand picture a boundary layer of 20 points. 



At the lateral and top boundary we introduce an additional bound- 
ary layer of nd points and w decreases from 1 to in this layer. We 
investigate different profiles regarding the weighting function, e.g. lin- 
ear, quadratic, cos and tanh (see figure 1). We investigate how the size 
of the boundary layer influences the solution. Table I and figure 3 show 
the result of our investigations. 
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Table I. Details of runs to reconstruct force-free and non 
force-free magnetic fields. All configurations have been cal- 
culated for a physical domain of 40 x 40 x 20 grid points. 
The first column specifies the profile of the boundary layer, 
the second column the size nd of the boundary layers, the 
third column the value of L in the computational box, 
the fourth column the value of Li in the physical domain 
and the fifth column the force-free condition (for force-free 
configuration) or the force balance (for MHS equilibria) av- 
eraged over the physical domain. We specify the start-error 
and discretisation error for each configuration. 
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4.1. Force-Free I 

Figure 3 top right panel shows that the field lines of a potential field 
reconstruction are clearly different from the original Low and Lou so- 
lution. This naturally leads to high values of the functional Lj and 
large J x B forces after the bottom boundary has been replaced by 
the original vector magnetogram. We first apply the optimization code 
without weighting function (nd = 0) . Here the boundaries of the physi- 
cal domain coincide with the computational boundaries. The lateral 
and top boundary have the value of the potential field during the 
iteration. Some low lying field lines are represented quite well (left 
hand picture in figure 3 second row). These field lines close to the box 
center are of course close to the bottom boundary and far away from 
the other boundaries. The (observed) bottom boundary has a higher 
influence on the field here than the potential lateral and top boundary. 
Other field lines, especially high reaching field lines deviate from the 
analytic solution (yellow dotted line). 

The values Lj and |J x B| provide a quantitative measure of the 
quality of the reconstructed magnetic field in the physical domain. High 
values correspond to a significant deviation from the force-free state. 
We applied the method of boundary relaxation (marked with b.r. in 
table I), but the result only slightly improved 2 . 

We investigate how the size and shape of a boundary layer influences 
the quality of the reconstruction. Both the comparison of the field lines 
(figure 3 panel 4 to 6) as well as the quantitative values in table I 
show that the quality of the reconstruction improves significantly with 
the size of the boundary layer (thickness in number of grid points nd). 
The larger computational box displaces the lateral and top boundary 
further away from the physical domain and consequently its influence 
on the solution decreases. As a result the magnetic field in the physical 
domain is dominated by the vector magnetogram data, which is exactly 
what is required. We find that a cos-profile of the weighting function 
provides the best results, followed closely by a tanh-profile. The main 
advantage of these profiles seem to be that they have smooth gradi- 
ents at the boundary of the physical domain to the boundary layer as 
well as at the boundaries of the computational box. We tried also a 
homogenous profile (w = 1 in the whole computational box, marked 
with horn in table I). In the homogeneous case w is equal 1 in the 
whole computational box and a boundary layer does not exist. The use 
of nd = 10 for the homogeneous case in table I is only for diagnostic 

2 The boundary relaxation method uses the iterative improvement — ^G for 
fields on the lateral and top boundary with G as in (16) in addition to (17) with 
w — 1. See ((Wiegelmann and Neukirch, 2003) section 4.2.3 for details. 
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Table II. The influence of noise. All con- 
figurations have been calculated for a 
physical domain of 40 x 40 x 20 grid points, 
a boundary layer of nd = 10 grid points 
and a cos profile in w. The first column 
specifies the noise level, the second col- 
umn the value of L in the computational 
box, the third column the value of Li in 
the physical domain and the forth column 
the force-free condition averaged over the 
physical domain. 



noise level 


L 


Li 


JxB 


[T 2 m] 


[T 2 m] 


[nNm- 3 ] 


no noise 


334 


65 


55 


1% 


678 


385 


150 


2% 
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1294 


267 


3% 


2847 


2440 
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5% 


6521 


6048 


505 


10% 


24240 


23332 


950 


15% 


52827 


51382 


1370 


20% 


92535 


90360 


1793 



reasons and indicates that Li and |J x B| have been computed in the 
same interior box (physical domain) as in the cases with weighting 
function for comparison. The effect that the lateral and top boundary 
are far away from the physical domain remains valid here, but the use 
of a weighting function in the boundary layers provides much better 
results. 

4.1.1. The influence of noise 

The previous calculations have been carried out under the assumption 
that the magnetic field on the boundary of the computational box is 
known exactly. Such an idealized situation will not be found when real 
vector magnetogram data are used. To keep control over the amount 
of uncertainty, we have carried out test runs by adding random noise 
to the vector magnetogram. We add the noise by multiplying the exact 
boundary conditions with a number 1 + 5 where 5 is a random number 
in the range — n\ < 8 < n\ and n% is the noise level. We investigate the 
effect of noise for Force-Free I with a boundary layer of nd = 10 grid 
points and a cos profile in w for different noise levels. 

Table II and figure 4 show our results. \/Z, \fL~l and |J x B| in- 
crease linearly with the noise level. The field line pictures in the upper 



opti04.tex; 4/02/2008; 1:13; p. 14 



Optimization code 



15 




Noiselevel Noiseleve 

Figure 4- Top row: The left hand side shows some field lines (same footpoints as in 
figure 3) for a noise level of 5% and the right hand side for a noise level of 10%. The 
yellow dotted field lines correspond to the original Low and Lou solution. We used 
a boundary layer of nd — 10 grid points and a cos profile in w. Bottom row: The 
left hand side shows a plot of the noise level against vLi an d the right hand side a 
plot of the noise level against |J x B|. The stars correspond to the data points and 
the line is a linear fit. 



panel of figure 4 show that low lying field lines are represented almost 
correctly while there are some deviations from the analytic Low and 
Lou solution for high reaching field lines. As the noise is completely 
random and independent between neighbouring grid points there are 
obvious difficulties regrading the discretisation. Our code uses forth 
order finite differences and consequently five grid points are required to 
compute derivatives. The highly oscillatory noise (The spatial variation 
of the noise corresponds to the spatial resolution of the grid.) naturally 
results in finite gradients which are linearly dependent from the noise 
level. | J x B| is linearly dependent from spatial deviations and does 
consequently also depend linearly on the noise level. The computation 
of L and Li contain gradients squared and are consequently quadratic 
dependent from the noise level. A pre-processing of raw magnetogram 
data, e.g. a Fourier filter or some smoothing might help to reduce the 
effect of random noise. 
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As the weighting function is designed to diminish the effect of the 
lateral and top boundary towards the solution, we investigate here how 
changes on these boundaries influence the solution. We undertook (for 
nd = 10 and a cos-profile) two runs with linear force free magnetic fields 
(aL x = 2.0 and aL x = —2.0) on the lateral and top boundaries. L,Lj 
and | J x B| are of the same order as for potential boundary conditions: 

potential : L = 334, U = 65, | J x B| = 55 

aL x = 2: L = 458, U = 92, | J x B| = 75 

aL x = -2 : L = 492, U = 62, | J x B| = 59 
The magnetic field lines calculated from these fields look exactly the 
same as for potential boundary conditions. We conclude that the influ- 
ence of the lateral and top boundary conditions towards the solution 
in the physical domain is indeed very small. 

4.2. Force-Free II 

The magnetogram on the right hand side of figure 2 is obviously bad 
conditioned. Significant parts of the magnetic flux are outside of the 
framed physical domain. In principle it would be better to always 
choose a magnetogram such that the majority of the flux is centered 
and the overall flux is approximately balanced but due to the limited 
size of vector magnetograms sometimes vector magnetograph data are 
not available over an entire active region. Here we investigate how this 
influences the quality of the reconstruction. We find that the recon- 
struction without weighting function provides much worse results in 
such a case than for a well centered force-free configuration. This is 
obvious caused by the potential field assumed on the lateral boundary 
in a region of a high magnetic flux and current density. The influence of 
the inconsistent lateral boundary conditions on the magnetic field in the 
physical domain is too strong. Here the method of boundary relaxation 
improves the quality of Li nearly by a factor of two, but the result 
still remains unsatisfactory. For runs with weighting function the result 
is basically similar as for Force-Free I. A large boundary layer (now 
including all active parts of the magnetogram) improves the quality of 
the reconstruction significantly. For the largest boundary layer nd = 20 
points wide the quality of the reconstruction is approximately equal to 
a 10 point wide boundary layer in the Force-Free I case. 

If for observational data only parts of an active region are available 
as vector magnetogram data a nonlinear reconstruction of the coronal 
magnetic field might become difficult or impossible. One could try to 
get the corresponding normal component of the photospheric magnetic 
field from other sources, e.g. the line of sight magnetograph MDI on 
SOHO and make assumptions regarding the transversal magnetic field. 
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Figure 5. Here we illustrate an example of a non force-free reconstruction. The 
left hand side shows a potential field reconstruction and the right hand side a 
MHS reconstruction with help of our optimization-code and a boundary layer of 
10 points with cos-profile. The yellow dotted lines show corresponding field lines for 
the exact solution. The colour-coding shows the normal magnetic field strength on 
the photosphere. 

The reconstructed magnetic field will then of course be influenced by 
these (not observed) assumptions. 

5. Tests for non force-free configurations 

For finite j3 configurations it is necessary to include pressure and gravity 
forces and solve the magneto hydro static equations (4)-(6). To test 
our code we use an analytic MHS equilibrium (see (Wiegelmann and 
Inhester, 2003) section 3.1.) with an average plasma (3 = 0.2. From this 
analytic solution we extract a photospheric vector magnetogram and 
the coronal density distribution. It is convenient to use the analytic 
density distribution here to test our code, but one has to keep in mind 
that for observational data the coronal density structure has to be 
reconstructed with help of a tomographic inversion. Corresponding ob- 
servational data (line of sight integrals of the coronal density structure) 
are expected from the Stereo-mission. The tomographic inversion of 
these data is a challenging problem on its own and we are not going 
to discuss it here. (Wiegelmann and Inhester, 2003) give an overview 
over the tomographic inversion algorithm. Here we consider the vector 
magnetogram and the plasma density distribution as given. The left 
hand panel of figure 5 shows a potential magnetic field reconstruction 
consistent with the normal component of the magnetic field extracted 
from an analytic magneto hydro static equilibrium. The comparison of 
the magnetic field lines of the potential field with the exact field show 
significant deviations. In table I we diagnose L and Lj similarly as for 
the force-free cases and we diagnose the force balance | J x B — VP| 
here. Similar as in the force-free case we find that the quality of the 
reconstruction is poor for an optimization without weighting function. 
A boundary layer of nd = 10 points together with a cos-profile of w 
improves the value of functional Li by a factor of 10 -2 in the physical 
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domain and the force balance 3 by a factor of about 10 . In the right 
hand panel of figure 5 we compare the reconstructed field with the 
analytic solution and both coincide within the graphical resolution. 
A larger boundary layer of nd = 20 further improves the quantita- 
tive measures of the non force-free reconstruction. The influence of 
a weighting function acts similarly for force-free and non force-free 
configuration. In both cases the reconstruction result is significantly 
improved. 

6. Conclusions 

In this paper we improved the optimization method for the recon- 
struction of non linear force-free and not force-free coronal magnetic 
fields. The optimization method minimizes a functional which consists 
of a quadratic form of the force balance and the solenoidal condition. 
While in previous optimization attempts the magnetic field needed to 
be described on all six boundaries of a computational box, our approach 
allows us to reconstruct the coronal magnetic field from the bottom 
boundary data alone. This is possible by the introduction of a boundary 
layer around the physical domain. The physical domain is a cubic area 
within which we want to reconstruct the coronal magnetic field consis- 
tent with photospheric vector magnetogram data. The boundary layer 
replaces the hard lateral and top boundary used previously. We showed 
that the limit of an infinitesimally thin boundary layer formally coin- 
cides with the original hard boundary. However, our test calculations 
show that a finite size weighted boundary yield much better results. 

We introduced a weighting function. In the physical domain the 
weighting function is unity. It drops monotonically in the boundary 
layer and reaches zero at the boundary of the computational box. At 
the boundary of the computational box we set the field to the value 
of the potential field computed from B n at the bottom boundary This 
choice is not only convenient due to its mathematic simplicity, but as 
well physical motivated. The coronal magnetic field is approaching a 
potential field in quite regions and high in the corona. The current 
program uses a cartesian geometry and is designed to reconstruct ba- 
sically isolated active regions. The photospheric magnetic field outside 
this active region (and a surrounding area) is ignored. This approach 
is the better justified the more the active region is isolated. Often, 
however, active regions are not always completely isolated but mag- 
netically connected with other active regions. To include this effect 

3 The functional Li contains the square of the force balance, see the definition 
11. 
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it would in principle be preferable to reconstruct the complete coronal 
magnetic field with vector magnetogram data on the whole photosphere 
as boundary. This would avoid the problems of prescribing the lateral 
boundaries and only the top boundary has to be chosen accordingly, 
similar to the source surface of potential field reconstructions. We do 
not expect any difficulties to apply our method to problems in spherical 
geometry. Unfortunately the required photospheric boundary data are 
not available because current vector magnetographs do not observe the 
entire solar surface. The situation will probably significantly improved 
by the vector spectromagnetograph of the the SOLIS project. It will 
deliver full disk vector magnetograms. 
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Appendix 
A. Derivation of F and G in (12). 

L= [ w(x,y,z) B 2 (n 2 a + n 2 b )d 3 x (A.l) 
Jv 

Cl a = B~ 2 [(V x B) x B + u] 
n h = B- 2 [(V • B) B] 

u = - M (VP + pW) (A.2) 

We vary L with respect to an iteration parameter t and get 

+ J w Sl b ■ J^[(V ■ B) B] <Px 

w (tf a + n 2 )B. — d 3 x (A.3) 

Our aim is now to use vector identities and Gauss law in such way 

that all terms contain a product with @M-. This will allow us to provide 
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explicit evolution equations for B to minimize L. The third term has 
the correct form already. We expand the first and second term 



1 dL 

2 ~dt 



+ 
+ 
+ 



/ w ci a ■ 

Jv 

/ w ft a ■ 

Jv 

/ w fl h ■ 

Jv 

w fib • 



dB 

(V x B) x 



x B 
dB 



d 3 x 



dt 



V 



dB\ 



B 



(V-B 



at ) 

dB 



d 3 x 
d 3 x 



dt 
dB 



- J vW (nl + nl)B. — d 3 x (A.4) 

The fourth and fifth term have the correct form. We apply the vector 
identities a • (b x c) = b • (c x a) = c • (a x b) to the first and second 



term 



2 ~dt 



w (v x ^ ) • (B x O a ) d 3 x 

V V ot 



f dB 

+ / w — ■ (O a x (V x B) d 3 x 
Jv dt 

f dB 
+ / w(n h -B)V- — d 3 x 

Jv dt 

/dB 
r w [fl b (V ■ B)] ■ — d 3 x 

- I 

Jv 



+ fig) B • ^ d 3 x 



(A.5) 



Term two, four and five have the correct form. We apply (V x a) • b = 
a • (V x b) + V • (a x b) to term 1 and • a = -a • + V • (atp) to 
term 3 



1 dL 

2 ~dt 



-x 



dB 

w -Tr ■ [V x (O a x B)] d 3 x 

\ at 

dB 

w V ■ {ft a x B) x — d 3 x 



dt 



dB 



+ / w — ■ (O a x (V x B) d A x 
Jv dt 

dB 



j v w V(!! b ' B) • — S 



+ 



w V 



v 



(n b .B,f 



d 3 x 
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+ J v w [fi b (V • B)] • ^ d 3 x 

w(n 2 a + n 2 b )B- — d 3 x 



(A.6) 



Until here the derivation has been identical with the method with- 
out weighting function (Wiegelmann and Inhester, 2003). Now we get 
additional terms with respect to the weighting function. We apply 
■0V • a = V • (ipa) — a • X7ip to term 2 and 5 



1 dL 

2 ~dt 



-I 

+ X 

X 

-x 
X 



<9B 

w -7T- • fV x (n a x B)l ci 3 x 
v at 

aB 1 

V • u; (fi a x B) x — d 3 x 



<9B 

(n a x b) x — 

<9B 



at _ 
• Vw d 3 x 



4 / u; — • (O a x (V x B)) d 6 x 
v at 

/9B 

w V(fi b • B) • — d 3 x 

ot 

aB 

io (fi b • B) — d 3 x 



+ / V- 



v 



( n b .B,f 



at 
<9B 



^ «, [n b (v-b)]. — d : 



(A.7) 



The terms 1,4,5,7,8 and 9 have the correct form. We apply Gauss' law 
to term 2 and 6 

I dL f dB 



2 dt 



+ X 



at 



[wVx(O a x B)] cf x 



w (fi a x B) x 



aB 



aB 

(O a X B) X — 



dt 
Vw d 3 x 



f r/B 

+ / — • ( w n a x (V x B)) d 3 x 
Jv at 

-x 

+ / fl( »n b .B).f # 



(9B 

«; V(ft b ■ B) • — d 3 x 
i at 

aB 
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<9B 



- J v [(fi b • B) Vw] ■ ^d 3 x 

r w [fi b (V • B)] • d 3 x 

- J v w(nl + nl)B.^d 3 x (A.8) 



Now all terms except the terms 2 and 3 have the correct form. We 
apply a • (b x c) = c • (a x b) to term 2 and 3 

1 dL f dB r /n xl 

= -|-.[.Vx( fia xB)]^ 



(9B 

[ u> n x (fi a x B)] • — d 2 x 
5B , 2 

x 



-Is 

f c)B 
+ / [Vw x (O a x B)] • — d 

jy at 

I' i9B 

+ / — • (O a x (V x B)) d 3 x 
Jv ot 

f (9B 

-X-v(n b -B).- gF A 

+ j n( u>n b B) —d 2 x 

/(9B 
«; [fib (V ■ B)] • — d 3 x 



Now all terms have the correct form and we write them more compact 

1 dL f 8B ~ - f 8B - , 2 

=>- — = -/ — — • F (I i - / — — • G d a; (A.10 

2 dt 7y at is 3* v 7 



F = wF + (O a xB)xVw + (f] b -B) Vw (A.ll) 

G = w G (A.12) 

F = V x (fi a xB)-!) a x(VxB) 

+V(fi b ■ B) - ft b (V ■ B) + {ill + «b) B ( A - 13 ) 

G = n x (n a x B) - n(ft b ■ B) (A.14) 
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